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Abstract 

We study a two-dimensional fluid of particles interacting through a spherically-symmetric and 
marginally soft two-body repulsion. This model can exist in three different crystal phases, one of 
them with square symmetry and the other two triangular. We show that, while the triangular solids 
first melt into a hexatic fluid, the square solid is directly transformed on heating into an isotropic 
fluid through a first-order transition, with no intermediate tetratic phase. In the low-pressure 
triangular and square crystals melting is reentrant provided the temperature is not too low, but 
without the necessity of two competing nearest-neighbor distances over a range of pressures. A 
whole spectrum of water-like fluid anomalies completes the picture for this model potential. 
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I. INTRODUCTION 



When confined to two-dimensional (2D) space, condensed matter behaves differently than 
in three dimensions (3D). A striking example is the physics of the Kosterlitz-Thouless tran- 
sition in superfiuid films, not to mention the phenomenon of high-temperature superconduc- 
tivity in the cuprates, which is intimately related to the motion of electrons within weakly- 
coupled copper-oxide layers. Another example is provided by 2D crystals, where thermal 
fluctuations are so strong as to rule out long-range translational (but not orientational) or- 
der for non-zero temperatures, leaving it open for the possibility of unconventional melting 
scenarios. The celebrated KTHNY theory of 2D melting thereby predicts a continuous two- 
stage melting from a crystalline to a hexatic phase and, subsequently, from the hexatic to an 
isotropic liquid [1]. The intermediate hexatic phase has short-range translational order (i.e., 
it is fluid-like) but quasi- long-range orientational order, characterized by a power-law decay 
of bond-angle pair correlations. Other possibilities are a standard first-order 2D-melting 
transition (as in Chui's theory [2]) or a stronger, i.e., first-order hexatic-to-isotropic tran- 
sition, as shown by Bernard and Krauth to be the case for hard discs [3j. In some cases, 
melting is one-stage and first-order but the hexatic phase can be accessed via a metastable 
route IHIS]. 

The KTHNY theory has been confirmed many times in simulation and experiment, espe- 
cially for particles with long-ranged interactions, though in a few cases with some controversy 
about the order of the two transition steps [5HT7] . From a computational point of view, the 
assessment of the order of 2D melting in practical cases can be hard, due to important 
finite-size effects and slow relaxation to equilibrium. An extreme case is hard discs where 
the exact nature of the melting transition could be established only recently, by employ- 
ing unprecedentedly large system sizes of the order of a million particles in a box of fixed 
volume |3j. Generally speaking, working in the isothermal-isobaric ensemble (rather than 
at constant volume) can ease the interpretation of the simulation data, since one cannot 
incur in any of the finite-size artifacts that harm constant-volume simulations. Recently, we 
have provided unambiguous evidence of the occurrence of continuous melting via a hexatic 
phase in the 2D Gaussian-core model (GCM) |T8]. This system, taken as a prototype of 
the phase behavior of "empty" polymers [I9_J, has long been known to exhibit reentrant 
melting (i.e., solid melting under isothermal compression) as well as water-like anomalies in 
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three dimensions j20]423] . In two dimensions, the GCM melting becomes two-staged, with 
an extremely narrow hexatic region whose properties comply with the predictions of the 
KTHNY theory. In particular, besides a reentrant isotropic fluid, also a reentrant hexatic 
phase exists. Recently, water-like anomalies have been reported even in the one-dimensional 
GCM [21], where however no properly deflned thermodynamic transition occurs. 

A class of systems which exhibit reentrant melting, solid polymorphism, and other water- 
like anomalies is formed by particles interacting via spherically-symmetric core-softened (CS) 
potentials For these systems, the strength of the repulsive component of the interac- 
tion undergoes some sort of weakening over a range of interparticle distances. The effects of 
particle-core softening on thermodynamic behavior were originally investigated by Hemmer 
and Stell [26] , who were interested in the possibility of multiple critical points and isostruc- 
tural solid-solid transitions in simple-fluid systems. A few years later. Young and Alder [27] 
showed that the phase diagram of the hard-core plus square-shoulder repulsion exhibits an 
anomalous melting line of the kind observed in Cs and Ce. Later, Debenedetti [2H] showed 
that systems of particles interacting via continuous CS potentials are capable of contracting 
when heated isobarically (a behavior which goes under the name of "density anomaly"). 
In the last decade or so, there has been a renewal of interest in the phase behavior of 
CS potentials, which has led to the discovery of many unusual properties such as multiple 
reentrant-melting lines, polymorphism both in the fluid and solid phases, stable cluster solids 
(for bounded interactions only), and a plethora of thermodynamic and dynamic water-like 
anomalies [25H1T] . A common feature of CS potentials is their ability to generate two dis- 
tinct length scales in the system, one related to the inner core and the other associated 
with the milder component of the repulsion [25]. Due to this, CS fluids are characterized 
by two competing, expanded and compact, local particle arrangements. Such an interplay 
has a disruptive influence on crystal order, the fluid phase being thus recovered upon com- 
pression. This two-scale or two-fluid property mimics the behavior of the more complex 
network-forming fluids, where the looser and denser local structures arise from the incessant 
building and breaking of the dynamic network generated by directional bonds. As a result, 
a paradigm has been established whereby the existence of two competing local structures is 
essential for the occurrence of anomalous behaviors in simple fluids. 

Unusual phase behavior has also been found in 2D systems. Scala and coworkers ^2] car- 
ried out molecular-dynamics simulations of the 2D square-shoulder plus square-well (SSSW) 
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potential finding water-like anomalies and identifying two different solid phases: a trian- 
gular crystal (low-density phase) and a square crystal (high-density phase). However, the 
observation of a hysteresis loop in their simulations suggested that the liquid-solid transition 
is actually first-order, thus washing out the possibility of a continuous melting to a hexatic 
phase. Wilding and Magee also studied the 2D SSSW potential and showed that the ther- 
modynamic anomalies of the model, rather than stemming from a metastable liquid-liquid 
critical point as previously surmised, were induced by the quasicontinuous nature of the 
2D freezing transition [43j. Almudallal and coworkers [H] contested this result, showing 
by various free-energy techniques that all transitions in the 2D SSSW model are in fact 
first-order. Malescio and Pellicane studied a 2D system of particles interacting through a 
potential consisting of an impenetrable hard core plus a square shoulder. They found a 
variety of stripe phases whose formation was imputed to the existence of two characteristic 
length scales [15]. All such investigations suggest that the existence of two length scales 
would be a prerequisite for observing water-like anomalies also in 2D. 

In a series of recent papers the two-fluid picture has been overthrown and a novel minimal 
scenario for the occurrence of anomalies has been established |16] - H9] . More specifically, it 
has been shown that a weakly-softened isotropic pair repulsion, with a single characteristic 
length which becomes more loosely defined in a range of pressures, is able to give rise to 
an unusual phase behavior. Our purpose here is to verify whether anything similar occurs 
in 2D, where in addition an interesting interplay with hexatic order may take place. Far 
from being purely academic, the present study can be relevant for many soft materials. 
For instance, one monolayer of A^-isopropylacrylamide (NIPA) microgel spheres confined 
between two glass cover slips is an ideal system to study 2D melting because the effective 
interparticle potential is short-ranged and repulsive, with a temperature-tunable volume 
fraction [HI [121 El EOl [51]. Systems like this would be natural candidates for detecting, by 
video optical microscopy or light scattering, the kind of phenomena that we are going to 
illustrate below. 

The plan of the paper is the following. After introducing our 2D model and method in 
Section II, we sketch the model phase diagram in Section III, highlighting the unconventional 
structure of the fluid phase which, similarly to other instances of CS repulsion, shows a 
number of water-like anomalies without any interplay between two characteristic length 
scales. The melting mechanism is studied in more detail in Section IV, where we show the 
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existence of a hexatic phase. Section V is finally devoted to concluding remarks. 

II. MODEL AND METHOD 

We consider a purely repulsive pair potential in two dimensions, modelled through an 
exponential form which was first introduced, about four decades ago, by Yoshida and Ka- 
makura (YK) [52]: 



where e and a set the energy and length scales, respectively, and a > 0. The YK potential 
behaves as for small r, and falls off very rapidly for large r. The smaller a the higher the 
degree of softness of 'u(r), i.e., the flatter the repulsive "shoulder" around r = a. Technically 
speaking [28], u{r) can be regarded as soft only for a ^ 2.3, since only in this range there 
is an interval of distances where the local virial function —ru'{r) decreases when the inter- 
particle separation decreases. However, it is over the wider range a ^ 5.5 that the phase 
diagram is expected to show a reentrant-fluid region [171 [55] . We shall here focus on the 
case a = 3.3 (Fig. 1), which was already shown to possess a rich anomalous phase behav- 
ior in three dimensions which cannot simply be explained by the existence of two distinct 
nearest-neighbor (NN) distances in the dense fluid [17]. 

We now briefly describe the methodology of the present investigation. After identifying 
the relevant crystal phases by means of total-energy calculations at zero temperature (of 
the type illustrated in, e.g., Ref. [M]), we explore the phase diagram by Monte Carlo (MC) 
simulation in the NPT ensemble. As is common practice, we adopt periodic boundary 
conditions and employ cell linked lists in order to speed up the simulation. While systems 
of about 1000 particles are suitable for investigating bulk properties and for determining the 
approximate location of phase boundaries, we consider larger systems of about 6000 particles 
for the search and characterization of the hexatic phase at selected pressures. For the same P 
values, we check the order of the melting transition independently through thermodynamic 
integration combined with exact free-energy calculations. For the fluid phase, a dilute gas 
is used as a reference state, whose chemical potential is calculated by the Widom method, 
whereas a low-temperature crystal is chosen as the starting point of the simulation in the 
solid region of the phase diagram. At this initial state, the excess Helmholtz free energy of 




(2.1) 
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the system is computed by the Einstein-crystal method [55] . 

III. FLUID STRUCTURE AND WATER-LIKE ANOMALIES 

For a = 3.3, the repulsive shoulder at cr, which is still visible in the plot of u{r) for 
a ~ 2, has by then completely faded out, surviving only in the form of a modest bump 
in the otherwise monotonously-decreasing profile of the virial function (see Fig. 1). Yet, 
this almost structureless potential shows three distinct stable crystal arrangements at T = 
in two dimensions, as it follows from the calculation of the chemical potential /i as a 
function of the pressure P for all five Bravais lattices and for the honeycomb lattice. Upon 
increasing P, the sequence of phases is triangular-square-triangular, the square crystal being 
thermodynamically stable in the pressure range 2.27-4.36, in units of e/a^ (from now on, 
pressure and temperature will be given in reduced units). 

Assuming that no other crystal phases come into play at nonzero temperatures, we first 
sketched the overall phase diagram by the heat-until-it-melts (HUIM) method. In practice, 
for selected P values we run a chain of MC simulations of the solid stable at the given P, 
for increasing T values at regular intervals of 0.0005, until we observe a clear jump in both 
the particle- number density p = N/V and total energy per particle E/N, which signal the 
melting of the solid. The character of this transition will be addressed in the next section. 
The HUIM method simply overlooks the possibility of solid overheating, but this is usually 
a reasonable approximation if one is not interested in very precise estimates of the melting 
temperature Tm {a posteriori, the typical error implied by the HUIM method in determining 
the melting point of the present system was about 10%, and smaller for the square crystal). 
The outcome of this analysis is reported in Fig. 2 where the most revealing feature is the 
reentrance of the fluid phase as the system, already settled in the low-density triangular 
solid, is further compressed at not too low temperatures. Fluid reentrance occurs even twice 
in a smaller range of temperatures above 0.020. 

As we already said in the Introduction, until very recently the conventional wisdom on the 
origin of the reentrant-melting phenomenon in CS systems with an unbounded interparticle 
repulsion rested on a competition, with destabilizing effects on crystal ordering, between 
two ways of arranging particles close to each other, which gives rise to either an expanded 
or a more compact structure in the dense fluid. This may only occur provided the pressure 
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is strong enough as to bring neighboring particles at distances close to a. Clearly, this 
explanation cannot work for a case like the present one, where it is hard to maintain that 
there are two different length scales in the potential. Rather, we may view the question from 
a solid-state perspective and argue that a succession of reentrant-melting lines in the phase 
diagram is simply the outcome of the existence of multiple solid phases at low temperature. 
On increasing pressure within the range of stability of any of such solids, the crystal strength 
first grows up to a maximum (and in parallel also T^), but then it progressively reduces 
on approaching the boundary of the next stable solid at higher pressure. After all, this is 
exactly the rationale behind the melting criterion discussed in Ref. [53], which in fact is very 
effective for soft repulsive interactions. 

We found confirmation that the present model shows only one repulsive length scale by 
computing the fluid radial distribution function (RDF) along the isotherm at T = 0.1, i.e., 
just above the maximum of the low-density triangular solid (see Fig. 3). Looking at 
the position of the main RDF peak, we see a systematic shift to lower and lower distances 
upon compression, with a slight widening of the peak around P ~ 2, i.e., next to the first 
reentrant-melting line. We interpret this as the evidence of a unique NN characteristic 
distance in the system, which becomes more loosely defined across the region of reentrant 
melting. This should be contrasted with what occurs for a more conventional CS repulsion, 
where the first RDF maximum is twin-peaked, with two definite NN distances that take 
turns at providing the absolute maximum for the RDF. 

Aside from the existence of one repulsive length scale rather than two, many water-like 
anomalies are also found in the present system, starting with a line of p maxima in the 
fluid close to the first reentrant-melting line (see Fig. 2). This is a typical occurrence in 
systems with CS potentials, where however the region of volumetric anomaly is usually 
more extended. We checked that the line of the density anomaly ceases to exist just before 
plunging into the solid region. Another type of anomaly is the so-called structural anomaly, 
that is a non-monotonous pressure behavior of the amount of "translational order" in the 
fluid, as measured through the value of minus the pair entropy per particle (— S2) [56] (see 
Fig. 2 and the left top panel of Fig. 4). Rather than monotonously increasing with pressure 
at constant temperature, the degree of spatial order reduces in the reentrant-fluid region, as 
a result of a looser definition of the NN distance. A non-sharp average separation between 
neighboring particles acts as a perturbing factor for the local order, bringing about a slight 
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decrease of — S2 with pressure. In the same range of pressures where —S2 decreases, the 
self-diffusion coefficient D, which we measure by a series of TVN molecular-dynamics runs, 
gets enhanced with pressure (Fig. 2 and right top panel of Fig. 4). Both lines of structural 
and diffusional anomalies appear to sprout out of the point of maximum for the low- 
density triangular solid and roughly terminate at the point of maximum Tm for the square 
solid. This is again similar to other CS systems, also for what concerns the crossing of 
the anomaly lines (see e.g. Fig. 11 of [57j). Finally, we checked for just one pressure value 
(P = 2) the existence of a temperature minimum in both the isobaric specific heat Cp and 
isothermal compressibility (lower panels of Fig. 4), similarly to what found for water at 
ambient pressure. Both minima fall in the high-temperature fluid region, much far from the 
maximum-density point for the same pressure. 

IV. HEXATIC BEHAVIOR 

The existence in our model of two solids with different crystal symmetry, i.e., triangular 
and square, gives the opportunity to investigate the intriguing possibility of two distinct 
intermediate phases (hexatic and tetratic) between the solid and the normal fluid, with the 
further bonus of a yet-to-be-observed hexatic-tetratic transition near the confluence point 
between the two solids and the fluid. 

In order to clarify the melting scenario for a given P, we run our computer code along 
two simulation paths, one starting from the solid phase at T = 0.005 and the other from a 
high-temperature fluid state. We advance in steps of AT = 0.005, equilibrating the system 
for long before generating an equilibrium trajectory of M sweeps (one sweep corresponding 
to trial MC moves), with M ranging between 5 x 10^ and 3 x 10^, depending on how 
far we are from melting. Once a guess of the transition point is made, we restart the 
simulation slightly ahead of it with a smaller AT and/or a larger M in order to better 
discriminate between first-order and continuous melting. Besides p and E/N, we measure 
two order parameters (OP), ipx and ipo, which are sensitive to the overall translational and 
orientational triangular/square order, respectively. The precise definition of both quantities 
has been given in Ref . [12] , with obvious modifications for the square-lattice case. Moreover, 
we keep track of the OP susceptibilities xt and xo, defined as N times the variance of 
the respective OP estimator. Finally, we calculate two orientational correlation functions 
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(OCF) [18], hQ{r) and h^ir), which inform on the typical size of a space region in which 
NN-bond angles are strongly correlated. The KTHNY theory predicts an algebraic r~^^^) 
large- distance decay of the OCF in the hexatic phase, at variance with what occurs in an 
isotropic fluid where the decay is much faster, i.e., exponential. According to the same 
theory, rj equals 1/4 at the transition point between hexatic and isotropic fluid. 

In Figs. 5-8, we report p and E/N for two different system sizes and various simulation 
protocols, as a function of T for P = 0.5, 2, 3, and 5. We clearly see that, while melting 
is continuous for P = 0.5, 2, and 5, it is certainly first-order for P = 3 (and 4, data not 
shown), as evidenced by the hysteresis loops. Based on our experience with the GCM [TS] . 
we can conjecture that there is a narrow region of hexatic phase for P = 0.5, 2, and 5 (to 
be confirmed later by the analysis of OPs and OCF), whereas no definite conclusion can be 
reached at this point for P = 3, where a tetratic phase can in principle exist even in presence 
of a first-order transition. A special remark is due for P = 5, where, in contrast with what 
occurs for smaller pressures, the crystal energy decreases upon heating, which means that 
the increase in kinetic energy is more than compensated for by the loss in potential energy, 
whose high rate of decrease is due to a NN distance lying in the harsh part of the potential 
core. 

For three pressures (0.5, 2, and 3), we plot the OPs and related susceptibilities in Figs. 9- 
11. For the first two pressures, we see that tpj- vanishes at a slightly smaller temperature 
than ipo = "ipe, which implies that the hexatic phase is confined to a narrow T interval not 
wider than 0.0005 for P = 0.5 (0.0015 for P = 2), as also witnessed by the maxima of 
the two susceptibilities occurring at slightly different T values. These temperature intervals 
compare well with the T range of the bridging region between the solid and fluid branches in 
Figs. 5 and 6. We thus confirm the same phenomenology of the GCM, namely that the width 
of the hexatic region increases with pressure. For P = 5, the findings are similar to P = 2, 
with roughly the same width of the hexatic region and comparable levels of orientational 
and translational order in the hot solid (data not shown). Going to P = 3, the picture is 
quite different since the two OPs now apparently vanish at the same temperature; at that 
point, the orientational susceptibility shows a spike (rather than a critical peak) which is 
usually associated with a first-order transition. All evidence suggests that standard first- 
order melting is a plausible explanation for P = 3, and the same conclusion can be made 
for P = 4 (data not shown). 
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A more direct evidence of the existence of a bond-angle ordered fluid (or a clue to its ab- 
sence) comes from the large-distance behavior of the OCF. We plot this function in Figs. 12- 
14 for various temperatures across the relevant region, for P — 0.5, 2, and 3, respectively. 
It appears that, for P = 0.5 and 2, h^^r) decays algebraically in a T region of limited ex- 
tent, roughly corresponding to the bridging region in Figs. 5 and 6. Moreover, the decay 
exponent in the hexatic region is smaller than 1/4, becoming larger only on transforming 
to the isotropic fluid, and the same is found for P = 5 (data not shown). On the contrary, 
for P = 3 and 4, /i4(r) switches directly from no decay at all to an exponential damping, 
showing that there is no tetratic phase in our model. 

The obvious question now arises as to why, at variance with the hexatic one, the tetratic 
phase is not stable in our model. This question is difficult to answer since it involves 
the consideration of the delicate equilibrium between the energy and entropy of the two 
competing phases, in this case the tetratic phase and the isotropic fluid phase. A possible 
hint could be the level of orientational order that is found in the two types of sohd shghtly 
before melting. If we look at Figs. 9-11, we see that the value of ■0o = 1/^4 ior P — 3 is less 
than a half of ipQ (note that the same holds for P = 4), and this would explain why, in the 
square-crystal case, long-range orientational order does not survive (even in the weakened 
form typical of a tetratic phase) the loss of quasi-long-range positional order determined by 
melting. 

Finally, we checked by an indipendent route the order of the melting transition for P = 
0.5, 2, 3, and 4. As anticipated, for each pressure we carried out MC simulations along two 
different paths, one beginning from a cold solid and the other from a low-density fluid. Using 
thermodynamic integration in combination with exact free-energy calculations at the initial 
points of the paths, we were able to obtain the system chemical potential /i along the solid 
and fluid branches as a function of T. For P = 0.5 and 2, we did not flnd any crossing of 
the two IJ>{T) curves, thus conflrming a continous melting transition. Notwithstanding the 
care we put in keeping under control any source of statistical error, we nevertheless found 
a small discrepancy (about 0.0005e, practically constant in a 0.01 wide T range across the 
melting transition) between the /j levels of the solid and the fluid, which is presumably due 
to the not-so-small AT employed along the paths far from melting. On the contrary, for 
P = 3 and 4, we observed a clear crossing between the two /x(T) curves, respectively at 
T — 0.0269 and T — 0.0176, which is consistent with the location of the density and energy 
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jumps (for P = 3, see Fig. 7), and suggestive of a first-order transition. 

Summing up, we collected multiple evidence of a continuous melting via a hexatic phase 
for the triangular crystals, while the melting transition is certainly discontinuous and "stan- 
dard" for the square crystal. Nothwithstanding the "small" sizes of the investigated samples, 
we think that our conclusions are robust since they result from many independent indicators 
of the phase-transition order. 

V. CONCLUSIONS 

We have analyzed the phase behavior of purely-repulsive 2D particles with a weakly- 
softened core. The nature of this repulsion is such as to determine a characteristic nearest- 
neighbor distance in the fluid phase whose statistical precision, expressed by the width 
of the main peak of the radial distribution function, shows a non-monotonous behavior 
with pressure at not too high temperatures. Notwithstanding the fact that the two-fluid 
paradigm of core-softened (CS) potentials does not apply here, we anyway observe the same 
phenomenology as in conventional CS systems, with solid polymorphism, multiple reentrant- 
melting lines, and many other water-like anomalies. While this is similar to the 3D case [17] , 
the melting transition of the present system is different, since it is continuous and occurs 
via an (even reentrant) hexatic phase for the triangular solids while being standard first- 
order for the intermediate square solid (i.e., no parallel tetratic phase exists). The hexatic 
behavior appears to be consistent with the KTHNY theory, as witnessed by the value of 
the decay exponent of the orientational correlation function at the hexatic-to-isotropic fluid 
transition temperature. 

Our findings could be relevant for many real soft-matter systems. Already nowa- 
days, monodisperse colloidal suspensions can be engineered in such a way as to exhibit 
a temperature-modulated repulsion with some amount of softness. Probably, the day is not 
far when by appropriate functionalization of the particle surface a colloidal system will be 
shown to exhibit the kind of phase-transition phenomena that we have illustrated here for 



11 



a rather specific model potential. 
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FIG. 1: (Color online). YK potential with a = 3.3 (Eq. (2.1), black solid line), its negative first- 
order derivative —u'{r) (dotted blue line), and the virial function —ru'{r) (dashed red line). 
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FIG. 2: (Color online). Phase diagram of the YK potential with a = 3.3. Two kind of melting 
points are shown: those determined by the HUIM method (blue open dots) and others (red-filled 
blue dots) marking either the onset of the hexatic phase {P = 0.5,2, and 5) or the crossing in 
temperature of solid and liquid chemical potentials (P = 3 and 4), as discussed in Section IV. 
Statistical errors are smaller than the symbols size. Another red dot is placed where the chemical 
potential of the low-density triangular (T) solid takes over that of the square (S) solid at T = 0.005. 
The black solid lines are schematic transition lines. The hexatic-isotropic fluid transition lines are 
not shown. The extent of the low-pressure hexatic region can be appreciated in the inset, which 
shows a magnification of the P interval from 1 to 2. Of similar width is the high-pressure hexatic 
region near P = 5. Further open symbols joined by straight lines mark the boundaries of anomaly 
regions: isothermal —S2 maxima and minima (left and right green squares); isothermal D minima 
and maxima (left and right magenta triangles); isobaric p maxima (red squares). We checked 
that the structural and diffusional anomalies are no longer present for T = 0.2 and T = 0.18, 
respectively; hence, we draw the dashed lines to mean that the two branches of each anomaly are 
actually sewed in from the top. 
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FIG. 3: YK potential with a = 3.3: radial distribution function g{r) for various pressures and for 
T = 0.08. Successive lines correspond to P = 0.2, 0.8, 1.6, 2.4, 3.2, 4. The arrow marks the direction 
along which the pressure increases. 




FIG. 4: Anomalous properties of the YK potential with a = 3.3. Left top panel: translational- 
order parameter —s2 (in units of ks) as a function of pressure for T = 0.08. Right top panel: 
self-diffusion coefficient D (units of (T(e/m)^/^, where m is the particle mass) for T = 0.08. Left 
bottom panel: isobaric specific heat Cp (units of fc^) as a function of temperature for P = 2. Right 
bottom panel: isothermal compressibility Kt (units of cr'^ /e) for P = 2. A clear minimum is seen 
in the plot of both response functions. 
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FIG. 5: (Color online). YK potential with a = 3.3 for P = 0.5: particle-number density (top) 
and total energy per particle (bottom) for two different sizes (A^ = 2688, red, green, and blue; 
N = 6048, black). Different colors denote different simulation protocols (red: AT = 0.005 and 
M = 5 X 10^; green: AT = 0.001 and M = 5 x 10^; blue: AT = 0.001 and M = 2 x 10^; black: 
AT = 0.0002 and M = 3 x 10^). Open dots and triangles refer to a heating path and a cooling 
path, respectively. 
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FIG. 6: (Color online). YK potential with a = 3.3 for P = 2: particle-number density (top) 
and total energy per particle (bottom) for two different sizes (A^ = 2688, red, green, and blue; 
N = 6048, black and yellow-filled black). Different colors denote different simulation protocols 
(red: AT = 0.005 and M = 5 x 10^; green: AT = 0.001 and M = 5 x 10^; blue: AT = 0.001 
and M = 2 X 10*^; black: AT = 0.001 and M = 3 x 10^; yellow-filled black: AT = 0.0003 and 
M = 3 X 10^). Open dots and triangles refer to a heating path and a cooling path, respectively. 
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FIG. 7: (Color online). YK potential with a = 3.3 for P = 3: particle-number density (top) 
and total energy per particle (bottom) for two different sizes (A^ = 2704, red, green, and blue; 
N = 6084, black). Different colors denote different simulation protocols (red: AT = 0.005 and 
M = 5 X 10^; green: AT = 0.001 and M = 5 x 10^; blue: AT = 0.001 and M = 2 x 10^; black: 
AT = 0.0001 and M = 3 x 10^). Open dots and triangles refer to a heating path and a cooling 
path, respectively. 
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FIG. 8: (Color online). YK potential with a = 3.3 for P = 5: particle-number density (top) 
and total energy per particle (bottom) for two different sizes (A^ = 2688, red, green, and blue; 
N = 6048, black). Different colors denote different simulation protocols (red: AT = 0.005 and 
M = 5 X 10^; green: AT = 0.001 and M = 5 x 10^; blue: AT = 0.001 and M = 2 x 10^; black: 
AT = 0.0005 and M = 2 x 10^). Open dots and triangles refer to a heating path and a cooling 
path, respectively. Inset, total energy per particle on a wider pressure range. 
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FIG. 9: YK potential with a = 3.3 for P = 0.5: order parameters and susceptibilities in the T 
range across the melting transition. Upper panels: the oricntational order parameter ipQ and its 
susceptibility xe- Dots and triangles mark data obtained by heating and by cooling, respectively. 
Lower panels: the translational order parameter tpx and its susceptibility xt on heating. All data 
from different protocols and sizes are reported, always preferring the most accurate estimate when 
more than one is available. 
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FIG. 10: YK potential with a = 3.3 for P = 2: order parameters and susceptibilities in the T 
range across the melting transition. See the caption of Fig. 9 for notation. 
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FIG. 11: YK potential with a = 3.3 for P = 3: order parameters and susceptibilities in the T 
range across the melting transition. See the caption of Fig. 9 for notation. 
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FIG. 12: (Color online). Orientational correlation function h^{r^ at selected temperatures across 
the hexatic region for P = 0.5 {N = 6048). Top: log-log plot; bottom: log-lin plot. Upon 
increasing T from 0.0598 to 0.0606 there is a qualitative change in the large-distance behavior of 
hQ^r), from constant (triangular solid) to power- low decay (hexatic fluid), up to exponential decay 
(isotropic fluid). Note that, consistently with the KTHNY theory, the decay exponent r] is less 
than 1/4 (which is the slope of the dashed straight line) in the hexatic phase. The slight recovery 
of correlations which is observed near the largest distance at which the OCF is computed (roughly 
corresponding to half of the simulation-box length) is a finite-size effect due to the use of periodic 
boundary conditions. 
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FIG. 13: (Color online). Orientational correlation function h^ir^ at selected temperatures across 
the hexatic region for P = 2 {N = 6048). Top: log-log plot; bottom: log-lin plot. Upon increasing 
T from 0.0345 to 0.0360 there is a qualitative change in the large-distance behavior of hQ{r), from 
constant (triangular solid) to power- low decay (hexatic fluid), up to exponential decay (isotropic 
fluid). Moreover, the decay exponent t] is less than 1/4 (which is the slope of the dashed straight 
line) in the hexatic phase. 
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FIG. 14: (Color online). Orientational correlation function /i4(r) at selected temperatures across 
the hexatic region for P = 3 (N = 6084). Top: log- log plot; bottom: log-lin plot. At variance with 
the triangular-lattice case, we assist to an abrupt change of decay mode as T goes from 0.0262 to 
0.0263, from constant (square solid) directly to exponential (isotropic fluid). 
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